MODULE CMAXEE
!this module solves the value function and policy functions.
  USE GLOBVAR
  USE NRUTIL, only : NR_BIG2, NR_SMALL, NR_SMALL2, ARRAY2D, NR_EPS
  use nr, only : sort
  USE MINIMIZATION, only : BFGS
  use simplex, only : Nelder_Meade
  USE INTERPOL
  use comf_trans, only : real2ab_exp, ab2real_exp
  use probability, only : CDF_Normal
  use COMFUNCS
  !use MAXINC
  use debugprint
  
  use MPI  !For Fortran 90, working with Intel version of openmpi.
  IMPLICIT NONE

CONTAINS
  
  !This is called in the head slaves, including the master if the master is a head slave
  !distribute nodes to slaves, collect results and then broadcast results to all slaves.
  SUBROUTINE cmaxee_VALUEFUNCTION()
     real(8) :: rcx(3), ev
     integer(4) :: i, iHstart, iHend, ierr, itemp, itemp2, numsent, status(MPI_STATUS_SIZE)
     
     ! So let's first solve the dynamic program from Period N_T, backwards
     do git = N_T, N_0, -1  !period
        git_pos = git - N_0 + 1 
        
        if (gcL_debugprint) write(*,'(A,I3,A,I3)') 'Solving color_mpi=',color_mpi,', git=',git
        
        
        !===========================================================
        ! Solve DP for a subset of H space in the head slave as well
        !===========================================================
        iHstart = gI_VTasks_mpi%a(new_id_mpi+1,1)
        iHend = gI_VTasks_mpi%a(new_id_mpi+1,2)
        if (iHstart <= iHend) then
           call solve_V4Ht(iHstart, iHend) 
        endif
        
       
        !======= MPI_Allgatherv results among the head slave and slaves
        if (new_numprocs_mpi>1) then
            !=== sync value and policy functions.
          iHstart = min(iHstart, N_H) 
          call MPI_Allgatherv(dvV_EV%a(1,1,1,1,1,iHstart,git_pos), gI_vsize_mpi*gI_VTasks_mpi%a(new_id_mpi+1,3), &
                      &  MPI_DOUBLE_PRECISION, dvV_EV%a(1,1,1,1,1,1,git_pos), gI_vsize_mpi*gI_VTasks_mpi%a(:,3), &
                      &  gI_vsize_mpi*gI_VTasks_mpi%a(:,4),    &
                      &  MPI_DOUBLE_PRECISION, new_comm_mpi, ierr)
          do i = 0, Ipt
              call MPI_Allgatherv(dcI(i)%a(1,1,1,1,1,iHstart,git_pos), gI_vsize_mpi*gI_VTasks_mpi%a(new_id_mpi+1,3),  &
                      &  MPI_DOUBLE_PRECISION, dcI(i)%a(1,1,1,1,1,1,git_pos), gI_vsize_mpi*gI_VTasks_mpi%a(:,3),  &
                      &  gI_vsize_mpi*gI_VTasks_mpi%a(:,4),     &
                      &  MPI_DOUBLE_PRECISION, new_comm_mpi, ierr)
              call MPI_Allgatherv(depsilon(i)%a(1,1,1,1,1,iHstart,git_pos), gI_vsize_mpi*gI_VTasks_mpi%a(new_id_mpi+1,3), & 
                      &  MPI_DOUBLE_PRECISION, depsilon(i)%a(1,1,1,1,1,1,git_pos), gI_vsize_mpi*gI_VTasks_mpi%a(:,3),  &
                      &  gI_vsize_mpi*gI_VTasks_mpi%a(:,4),     &
                      &  MPI_DOUBLE_PRECISION, new_comm_mpi, ierr)
          enddo !do i
          do i = 0, Ipt+1
             call MPI_Allgatherv(dcc(i)%a(1,1,1,1,1,iHstart,git_pos), gI_vsize_mpi*gI_VTasks_mpi%a(new_id_mpi+1,3),  &
                      &  MPI_DOUBLE_PRECISION, dcc(i)%a(1,1,1,1,1,1,git_pos), gI_vsize_mpi*gI_VTasks_mpi%a(:,3), &
                      &  gI_vsize_mpi*gI_VTasks_mpi%a(:,4),     &
                      &  MPI_DOUBLE_PRECISION, new_comm_mpi, ierr)
             call MPI_Allgatherv(dvV(i)%a(1,1,1,1,1,iHstart,git_pos), gI_vsize_mpi*gI_VTasks_mpi%a(new_id_mpi+1,3), &
                      &  MPI_DOUBLE_PRECISION, dvV(i)%a(1,1,1,1,1,1,git_pos), gI_vsize_mpi*gI_VTasks_mpi%a(:,3), &
                      &  gI_vsize_mpi*gI_VTasks_mpi%a(:,4),     &
                      &  MPI_DOUBLE_PRECISION, new_comm_mpi, ierr)
             call MPI_Allgatherv(dciRECSS(i)%a(1,1,1,1,1,iHstart,git_pos), gI_vsize_mpi*gI_VTasks_mpi%a(new_id_mpi+1,3), &
                      &  MPI_INTEGER, dciRECSS(i)%a(1,1,1,1,1,1,git_pos), gI_vsize_mpi*gI_VTasks_mpi%a(:,3), &
                      &  gI_vsize_mpi*gI_VTasks_mpi%a(:,4),     &
                      &  MPI_INTEGER, new_comm_mpi, ierr)
          enddo !do i
        endif !if (new_numprocs_mpi>1) then
        

     enddo !git
        
     !===========================================================================
     ! print out value and policy functions for all grids. This is for debugging
     ! !!!! WARNING: the output file is very large. !!!!!!
     !===========================================================================
     !if (gcL_printvalue) call printvalues(color_mpi)

     if (gcL_debugprint) write(*,'(A,I4,A,I2)') "iter#", gcI_itercount, &
         & ', finished value/policy function solving in color of ',color_mpi
          
     RETURN
  END SUBROUTINE cmaxee_VALUEFUNCTION

  !=======================================================================================================
  !END OF calculate value and optimal consumption for retiree at period git (j=N_T-git)
  !=======================================================================================================
        
!=======================================================================================================
!solve value/policy function for grid giH at time git (slaves and head slaves)
!=======================================================================================================
subroutine solve_V4Ht(istart, iend)
  
 implicit none
 integer, intent(in) :: istart, iend
 integer :: i
 real(8) :: ev,rcx(3), rptx, rVy(5)
 
 if (istart>N_H .OR. istart>iend) then
     write(*,*) "EXCEED!!!" !for debug only
     return
 endif
  
 do giH = istart, iend, 1
  gH = sgridH%a(git_pos, giH) !human capital
  do giRECSS = sgridRECSS%d1,1,-1
     if ((giRECSS == 2 .and. git <= pERA) .or. (giRECSS == 1 .and. git>=71 .and. pNRA<71)) then
        !cannot receive ss benefit before pERA, so no state with ss=1 on or before pERA.
        !No DRC is given for retirement at age 70 or older, so after 70 everyone takes ss.
        !giRECSS represents if receiving SS benefit in the previous period.
        !giRECSS==1 and giRECSS_next==2 means receiving SS benefit the first time at Period t.
        do giFS = 1, N_FS  !family status
           
           do i = 0, Ipt+1 
               dvV(i)%a(:, :, :, giRECSS, giFS, giH, git_pos) = 0.0d0
           enddo !do i
           
           dvV_EV%a(:, :, :, giRECSS, giFS, giH, git_pos) = 0.0d0
           
        enddo !do giFS
        cycle
        !The cycle statement will move at once to the next iteration
        !without further progress through the loop body for the current iteration
     endif
     do giHEALTH = 1, N_HEALTH   !health status 1-excellent, 2-good, 3-bad, 4-disabled
        !IHasHealthAge>=2: has health status starting from age IHasHealthAge
        !when git<IHasHealthAge, giHEALTH should be always 2 (good). 
        if ( (IHasHealthAge>0) .AND. (git<IHasHealthAge) .AND. (giHEALTH .NE. 2) ) then
           do giFS = 1, N_FS  !family status
              
              do i = 0, Ipt+1 
                  dvV(i)%a(:, :, giHEALTH, giRECSS, giFS, giH, git_pos) = 0.0d0
              enddo !do i
              
              dvV_EV%a(:, :, giHEALTH, giRECSS, giFS, giH, git_pos) = 0.0d0
              
           enddo !do giFS
           cycle
        endif   !if ( (IHasHealthAge>=2)
        do giFS = 1, N_FS  !family status  
          do giAIME = sgridAIME%d2, 1, -1
             gAIME = sgridAIME%a(git_pos, giAIME) !AIME
             do giA = 1, VN_A(git)   !Asset. Must start from 1, the lower end (involving lower bound)
                gA = sgridA%a(git_pos,giA) ! Asset
                
                !======================the innerest loop====================================
                !choose (c, I, L, recss) to maximize value function
                !===========================================================================
                if (giRECSS == 2) then  !already receiving ss.
                   giRECSS_next = 2  !receiving ss is absorbing. the filter makes sure that now git >= pERA
                   !case with leisure = 0
                   call dcmaxee_getVAL(dvV(0)%a(giA, giAIME, giHEALTH, giRECSS, giFS, giH, git_pos), &
                                 &  dcc(0)%a(giA, giAIME, giHEALTH, giRECSS, giFS, giH, git_pos), &
                                 &  dcI(0)%a(giA, giAIME, giHEALTH, giRECSS, giFS, giH, git_pos), 0)
                   dciRECSS(0)%a(giA, giAIME, giHEALTH, giRECSS, giFS, giH, git_pos) = giRECSS_next-1
                   !case with leisure = 1
                   call dcmaxee_getVAL(dvV(1)%a(giA, giAIME, giHEALTH, giRECSS, giFS, giH, git_pos), &
                                 &  dcc(1)%a(giA, giAIME, giHEALTH, giRECSS, giFS, giH, git_pos), &
                                 &  rcx(1), 1)
                   dciRECSS(1)%a(giA, giAIME, giHEALTH, giRECSS, giFS, giH, git_pos) = giRECSS_next-1
                   if (Ipt == 1) then
                      !case with leisure = 2 (part time)
                      call dcmaxee_getVAL(dvV(2)%a(giA, giAIME, giHEALTH, giRECSS, giFS, giH, git_pos), &
                                 &  dcc(2)%a(giA, giAIME, giHEALTH, giRECSS, giFS, giH, git_pos), &
                                 &  dcI(1)%a(giA, giAIME, giHEALTH, giRECSS, giFS, giH, git_pos), 2)
                      dciRECSS(2)%a(giA, giAIME, giHEALTH, giRECSS, giFS, giH, git_pos) = giRECSS_next-1 
                   endif
                   
                   
                else ! not receiving ss yet
                   if (git<70) then
                      !===== case with leisure = 0
                      giRECSS_next = 1 !not receive ss at git
                      call dcmaxee_getVAL(dvV(0)%a(giA, giAIME, giHEALTH, giRECSS, giFS, giH, git_pos),  &
                                 &  dcc(0)%a(giA, giAIME, giHEALTH, giRECSS, giFS, giH, git_pos), &
                                 &  dcI(0)%a(giA, giAIME, giHEALTH, giRECSS, giFS, giH, git_pos), 0)
                      dciRECSS(0)%a(giA, giAIME, giHEALTH, giRECSS, giFS, giH, git_pos) = giRECSS_next-1
                      if (git >= pERA) then ! check if they should start collecting social security benefit at this period
                         giRECSS_next = 2
                         ev = - NR_BIG2
                         call dcmaxee_getVAL(ev, rcx(1), rcx(2), 0)
                         if (ev >= dvV(0)%a(giA, giAIME, giHEALTH, giRECSS, giFS, giH, git_pos)) then
                            dcc(0)%a(giA, giAIME, giHEALTH, giRECSS, giFS, giH, git_pos) = rcx(1)
                            dcI(0)%a(giA, giAIME, giHEALTH, giRECSS, giFS, giH, git_pos) = rcx(2)
                            dciRECSS(0)%a(giA, giAIME, giHEALTH, giRECSS, giFS, giH, git_pos) = giRECSS_next-1
                            dvV(0)%a(giA, giAIME, giHEALTH, giRECSS, giFS, giH, git_pos) = ev
                         endif
                      endif !if (git >= pERA)
                      !=====case with leisure = 1
                      giRECSS_next = 1 !not receive ss at git
                      call dcmaxee_getVAL(dvV(1)%a(giA, giAIME, giHEALTH, giRECSS, giFS, giH, git_pos),  &
                                    &  dcc(1)%a(giA, giAIME, giHEALTH, giRECSS, giFS, giH, git_pos), &
                                    &  rcx(1), 1)
                      dciRECSS(1)%a(giA, giAIME, giHEALTH, giRECSS, giFS, giH, git_pos) = giRECSS_next-1
                      if (git >= pERA) then
                         giRECSS_next = 2
                         ev = - NR_BIG2
                         call dcmaxee_getVAL(ev, rcx(1), rcx(2), 1)
                         if (ev >= dvV(1)%a(giA, giAIME, giHEALTH, giRECSS, giFS, giH, git_pos)) then
                            dcc(1)%a(giA, giAIME, giHEALTH, giRECSS, giFS, giH, git_pos) = rcx(1)
                            dciRECSS(1)%a(giA, giAIME, giHEALTH, giRECSS, giFS, giH, git_pos) = giRECSS_next-1
                            dvV(1)%a(giA, giAIME, giHEALTH, giRECSS, giFS, giH, git_pos) = ev
                         endif
                      endif !if (git >= pERA)
                      !===== case with leisure = 2 PT (part time)
                      if (Ipt == 1) then                         
                         giRECSS_next = 1 !not receive ss at git
                         call dcmaxee_getVAL(dvV(2)%a(giA, giAIME, giHEALTH, giRECSS, giFS, giH, git_pos),  &
                                 &  dcc(2)%a(giA, giAIME, giHEALTH, giRECSS, giFS, giH, git_pos), &
                                 &  dcI(1)%a(giA, giAIME, giHEALTH, giRECSS, giFS, giH, git_pos), 2)
                         dciRECSS(2)%a(giA, giAIME, giHEALTH, giRECSS, giFS, giH, git_pos) = giRECSS_next-1
                         if (git >= pERA) then ! check if they should start collecting social security benefit at this period
                            giRECSS_next = 2
                            ev = - NR_BIG2
                            call dcmaxee_getVAL(ev, rcx(1), rcx(2), 2)
                            if (ev >= dvV(2)%a(giA, giAIME, giHEALTH, giRECSS, giFS, giH, git_pos)) then
                               dcc(2)%a(giA, giAIME, giHEALTH, giRECSS, giFS, giH, git_pos) = rcx(1)
                               dcI(1)%a(giA, giAIME, giHEALTH, giRECSS, giFS, giH, git_pos) = rcx(2)
                               dciRECSS(2)%a(giA, giAIME, giHEALTH, giRECSS, giFS, giH, git_pos) = giRECSS_next-1
                               dvV(2)%a(giA, giAIME, giHEALTH, giRECSS, giFS, giH, git_pos) = ev
                            endif
                         endif !if (git >= pERA) 
                      endif   !if (Ipt == 1)
                   else  !git>=70. Always take ss
                      !===== case with leisure = 0
                      giRECSS_next = 2 !always take ss
                      call dcmaxee_getVAL(dvV(0)%a(giA, giAIME, giHEALTH, giRECSS, giFS, giH, git_pos),  &
                                     &  dcc(0)%a(giA, giAIME, giHEALTH, giRECSS, giFS, giH, git_pos), &
                                     &  dcI(0)%a(giA, giAIME, giHEALTH, giRECSS, giFS, giH, git_pos), 0)
                      dciRECSS(0)%a(giA, giAIME, giHEALTH, giRECSS, giFS, giH, git_pos) = giRECSS_next-1
                      !=====case with leisure = 1
                      giRECSS_next = 2 !always take ss
                      call dcmaxee_getVAL(dvV(1)%a(giA, giAIME, giHEALTH, giRECSS, giFS, giH, git_pos),  &
                               &  dcc(1)%a(giA, giAIME, giHEALTH, giRECSS, giFS, giH, git_pos), rcx(1), 1)
                      dciRECSS(1)%a(giA, giAIME, giHEALTH, giRECSS, giFS, giH, git_pos) = giRECSS_next-1
                      !===== case with leisure = 0
                      if (Ipt == 1) then
                         giRECSS_next = 2 !always take ss
                         call dcmaxee_getVAL(dvV(2)%a(giA, giAIME, giHEALTH, giRECSS, giFS, giH, git_pos),  &
                                     &  dcc(2)%a(giA, giAIME, giHEALTH, giRECSS, giFS, giH, git_pos), &
                                     &  dcI(1)%a(giA, giAIME, giHEALTH, giRECSS, giFS, giH, git_pos), 2)
                         dciRECSS(2)%a(giA, giAIME, giHEALTH, giRECSS, giFS, giH, git_pos) = giRECSS_next-1 
                      endif
                   endif !if (git<70)
                endif !if (giRECSS == 2)
                !========================End of the innerest loop=================================
                
             enddo  !giA
             
             !call cmaxee_consistency_VCinA()
             
             !calculate Emax and epsilon
             rVy(3) = gr_leisure_taste(git,giHEALTH,giFS)
                          
             do giA = 1, VN_A(git)
                rcx(1) = dvV(0)%a(giA, giAIME, giHEALTH, giRECSS, giFS, giH, git_pos) !leisure=0
                rcx(2) = dvV(1)%a(giA, giAIME, giHEALTH, giRECSS, giFS, giH, git_pos) !leisure=1
                if (Ipt == 0) then
                    if (abs(pa_epsilon) <= 0.0d0) then  !no stochastic shock
                       rcx(3) = exp(rVy(3))
                       rcx(2) = rcx(2) + rcx(3)
                       if (rcx(2) >= rcx(1)) then
                          rVy(1) = 0.0d0  !prob of (leisure=0)
                          depsilon(0)%a(giA, giAIME, giHEALTH, giRECSS, giFS, giH, git_pos) = - 1.0d3 !leisure=1 for sure
                       else
                          rVy(1) = 1.0d0  !prob of (leisure=0)
                          depsilon(0)%a(giA, giAIME, giHEALTH, giRECSS, giFS, giH, git_pos) = 1.0d3  !leisure=0 for sure
                       endif
                       dvV_EV%a(giA, giAIME, giHEALTH, giRECSS, giFS, giH, git_pos) = & 
                           &  rVy(1) * rcx(1) + (1.0d0-rVy(1)) * rcx(2) 
                    else   !if (pa_epsilon>0.0d0) then
                       rcx(3) = rcx(1) - rcx(2)  !leisure 0 - leisure 1
                       if (rcx(3) < NR_EPS) rcx(3) = NR_EPS
                       rcx(3) = (log(rcx(3))-rVy(3)-grUVmultiplier_ln)/pa_epsilon  !depsilon, with standard normal distribution
                       depsilon(0)%a(giA, giAIME, giHEALTH, giRECSS, giFS, giH, git_pos) = rcx(3)
                       rVy(1) = CDF_Normal(rcx(3), 0.0d0, 1.0d0) !Phi(depsilon)
                       
                       rcx(3) = gr_leisure_condexp(git,giHEALTH,giFS) * CDF_Normal(pa_epsilon - rcx(3), 0.0d0, 1.0d0) 
                                          
                       dvV_EV%a(giA, giAIME, giHEALTH, giRECSS, giFS, giH, git_pos) = & 
                           &  rVy(1) * rcx(1) + (1.0d0 - rVy(1)) * rcx(2) + rcx(3)
                    endif  !if (abs(pa_epsilon)<NR_EPS) then
                else   !Ipt == 1 !assuming there is always stochastic shock
                   rptx = dvV(2)%a(giA, giAIME, giHEALTH, giRECSS, giFS, giH, git_pos) !leisure==2 (part time)
                   
                   rcx(3) = rcx(1) - rcx(2)  !leisure 0 - leisure 1
                   rVy(1) = min((rcx(1)-rptx)/max(NR_EPS,plpt_t(git,giFS)), rcx(3))   !gamma_tmin. Condition A8. Leisure 0 - p
                   rVy(2) = max((rptx-rcx(2))/max(NR_EPS,(1.0d0-plpt_t(git,giFS))), rcx(3))  !gamma_tmax. Leisure p - 1
                   
                   
                   rVy(1) = (log(max(NR_EPS, rVy(1)))-rVy(3)-grUVmultiplier_ln)/pa_epsilon !epsilon_tmin. /pa_epsilon is to get a standard normal with sd=1
                   rVy(2) = (log(max(NR_EPS, rVy(2)))-rVy(3)-grUVmultiplier_ln)/pa_epsilon !epsilon_tmax
                   depsilon(1)%a(giA, giAIME, giHEALTH, giRECSS, giFS, giH, git_pos) = rVy(1)  !epsilon_tmin
                   depsilon(0)%a(giA, giAIME, giHEALTH, giRECSS, giFS, giH, git_pos) = rVy(2)   !epsilon_tmax
                   
                   rVy(4) = CDF_Normal(rVy(1), 0.0d0, 1.0d0)   !Phi(epsilon_tmin)
                   rVy(5) = CDF_Normal(rVy(2), 0.0d0, 1.0d0)   !Phi(epsilon_tmax)
                   
                   rcx(3) = rVy(4) * rcx(1) + (1.0d0 - rVy(5)) * rcx(2)  &
                          & + gr_leisure_condexp(git,giHEALTH,giFS) * CDF_Normal(pa_epsilon-rVy(2),0.0d0,1.0d0)
                   
                   if (rVy(4) < rVy(5)) then
                      rcx(3) = rcx(3) + (rVy(5)-rVy(4))*rptx   &
                          & + gr_leisure_condexp(git,giHEALTH,giFS) * ( CDF_Normal(pa_epsilon-rVy(1),0.0d0,1.0d0)  &
                          &   - CDF_Normal(pa_epsilon-rVy(2), 0.0d0, 1.0d0) )
                   endif !if (rVy(4) < rVy(5))
                   dvV_EV%a(giA, giAIME, giHEALTH, giRECSS, giFS, giH, git_pos) = rcx(3)
                endif  !if (Ipt == 0) 
             enddo !do giA = 1, 
             
             !It is important that at the lower end, dvV_EV is strictly increasing with A, to avoid the extrapolation problem.
             ! dvV_EV should be concave.
             call cmaxee_consistency_EVinA()
             
          enddo !giAIME
       enddo  !giFS
    enddo    !giHEALTH 
  enddo  !giRECSS
 enddo  !giH
 return
end subroutine solve_V4Ht
!=======================================================================================================
!END OF solve value/policy function for grid giH at time git
!=======================================================================================================

  !=======================================================================================================
  !calculate value and optimal consumption for period git
  !need gAIME, gA, gH, giHEALTH
  ! only for discrete leisure case
  !=======================================================================================================
  SUBROUTINE dcmaxee_getVAL(ev_, c_, I_, iL_)
    IMPLICIT NONE
    real(8), INTENT(INOUT) :: ev_
    integer, intent(IN) :: iL_
    real(8), INTENT(OUT) :: c_, I_
    
    real(8) :: rcx2(2), Vrcx1(1), ev, rtemp
    integer :: i,j, itemp
    
    gI_p_L = iL_
    if (iL_ == 2) then  !pt
        gr_p_L = 0.5d0
    else
        gr_p_L = dble(iL_)
    endif
    
    if (git==N_T) then
       ev_ = cmaxee_valfun_pT(iL_)
       c_ = gr_cc
       I_ = 0.0d0  !!**last period. I = 0, regardless of L
    else !git<N_T
       !upper bound of the C, git<N_T here
       rtemp = pw(git) * gH * (1.0d0-gr_p_L)  !labor income if working full hours (full time or part time)
              
       call A_AIME_prime( 0.0d0, rtemp, rcx2)
       gr_p_cbound(2) = (rcx2(1)-sgridA%a(git_pos+1,1))*pr_inv  !max of c such that asset next period is in the lower bound

       !check if the case where a positive amount of consumption is not affordable
       if (gr_p_cbound(2) <= NR_EPS) then  !cannot afford any consumption. sgridA does not depend on H
          c_ = grCF
          I_ = 0.0d0
          Vrcx1(1) = ab2real_exp(gr_p_Ibound(iL_,1), gr_p_Ibound(iL_,2), I_) !I = 0.0d0
          ev_ = - cmaxee_valfun_I(Vrcx1)
          return
       endif  !if (iL_==1 .AND. gr_p_cbound(2) <= 0.0d0)
       
       !lower bound of c, initial guess of c and I
       ev = gr_BL * 0.5d0
       if (giA == 1) then
          rcx2(1) = ev !initial c
          if ((iL_ .NE. 1) .AND. (gcI_robust<8 .OR. gcI_robust>13)) then
             rcx2(2) = min(0.1d0, 0.002d0 * dble(N_T-git))  !initial guess of I
          else
             rcx2(2) = 0.0d0  !leisure=1 or exogeneous models
          endif
       else
          rtemp = dcc(iL_)%a(giA-1,giAIME,giHEALTH,giRECSS,giFS,giH,git_pos) 
          rcx2(1) = rtemp  !initial guess of c
          
          if ((iL_ .NE. 1) .AND. (gcI_robust<8 .OR. gcI_robust>13)) then
              if (iL_==0) then  !FT
                  rcx2(2) = dcI(0)%a(giA-1,giAIME,giHEALTH,giRECSS,giFS,giH,git_pos)      !initial guess of I
              else  !PT
                  rcx2(2) = dcI(1)%a(giA-1,giAIME,giHEALTH,giRECSS,giFS,giH,git_pos)      !initial guess of I
              endif
              
          else  !leisure=1 or exogeneous models
             rcx2(2) = 0.0d0 
          endif
       endif       
              
       !fine tune the initial c and I
       if (gr_p_cbound(2)-gr_p_cbound(1) >= gr_BL) then
          rcx2(1) = min(rcx2(1), 50.0d0)  !initial c
          rcx2(1) = max(gr_p_cbound(1)+NR_EPS, min(gr_p_cbound(2)-NR_EPS, rcx2(1)))
       else
          rcx2(1) = (gr_p_cbound(1)+gr_p_cbound(2)) * 0.5d0 
       endif       
       rcx2(1) = ab2real_exp(gr_p_cbound(1), gr_p_cbound(2), rcx2(1))  !c
       
       if ((iL_ .NE. 1) .AND. (gcI_robust<8 .OR. gcI_robust>13)) then  !leisure = 0 or 2 in the normal BP model. Search I and c
          rcx2(2) = max(1.0d-5, min(0.3d0, rcx2(2)))   !initial guess of I
          rcx2(2) = ab2real_exp(gr_p_Ibound(iL_,1), gr_p_Ibound(iL_,2), rcx2(2))  !I
          
             !search over c and I simultaneously, using Simplex
                          
             call Nelder_Meade(rcx2, gcF_vtol, cmaxee_valfun_cI, 1, gcI_imax)
             
             ev_ = - cmaxee_valfun_cI(rcx2)       
             c_ = real2ab_exp(gr_p_cbound(1), gr_p_cbound(2), rcx2(1))
             I_ = real2ab_exp(gr_p_Ibound(iL_,1), gr_p_Ibound(iL_,2), rcx2(2)) 
             
             !check c=grCF case
             if ( min(gA,c_) < 20.0d0 ) then
                Vrcx1(1) = ab2real_exp(gr_p_Ibound(iL_,1), gr_p_Ibound(iL_,2), I_) !I
                call Nelder_Meade(Vrcx1, gcF_vtol, cmaxee_valfun_I, 1, gcI_imax)
                ev = - cmaxee_valfun_I(Vrcx1)
                if (ev > ev_) then
                   c_ = grCF
                   I_ = real2ab_exp(gr_p_Ibound(iL_,1), gr_p_Ibound(iL_,2), Vrcx1(1))
                   ev_ = ev
                endif
             endif          
          
       else   !leisure==1 or no need for I
          gr_p_I = 0.0d0     !since leisure = 1.0d0, or no need for I
          Vrcx1(1) = rcx2(1)  !initial c
          call Nelder_Meade(Vrcx1, gcF_vtol, cmaxee_valfun_c, 1, gcI_imax)
          
          ev_ = - cmaxee_valfun_c(Vrcx1)
          c_ = real2ab_exp(gr_p_cbound(1), gr_p_cbound(2), Vrcx1(1))
          I_ = 0.0d0
          
          !check c=grCF case
          if ( min(gA,c_) < 20.0d0 ) then
                Vrcx1(1) = ab2real_exp(gr_p_Ibound(iL_,1), gr_p_Ibound(iL_,2), I_) !I = 0.0d0
                ev = - cmaxee_valfun_I(Vrcx1)
                if (ev > ev_) then
                   c_ = grCF
                   I_ = 0.0d0
                   ev_ = ev
                endif
          endif  !if ( min(gA,c_) < 20.0d0 )
       endif  !if (iL_ == 0 .AND. (gcI_robust<8 .OR. gcI_robust>11))       
    endif   !if (git==N_T)
        
    return
    
  end subroutine dcmaxee_getVAL
 
  !=========================================================
  ! Calculate value function given c & I, and L=gr_p_L
  !=========================================================
  real(8) function cmaxee_valfun_cI(rcx_)
    !calculate the value function given c & I for Period git
    ! Since it is to be minimized, flip the sign
    REAL(8), INTENT(IN) :: rcx_(:)
    real(8) :: c_, I_, L_, rtemp, rx2(2)
    c_ = real2ab_exp(gr_p_cbound(1), gr_p_cbound(2), rcx_(1))
    I_ = real2ab_exp(gr_p_Ibound(gI_p_L,1), gr_p_Ibound(gI_p_L,2), rcx_(2))
    L_ = gr_p_L    
    IF ((I_ .NE. I_) .or. (c_ .NE. c_)) THEN	! this will be true if I_ or L_ is NaN
       cmaxee_valfun_cI = NR_BIG2
       return
    endif
    if ( (I_ < 0.0d0) .OR. (L_ < 0.0d0) .OR. ((I_+L_) > 1.0d0+NR_EPS) ) then
       cmaxee_valfun_cI = NR_BIG2
       return
    endif
    !The next period state variables
    if (gI_p_L==1) then !leisure = 1
       rtemp = 0.0d0
    elseif (gI_p_L==0) then  !leisure=0
       rtemp = pw(git) * gH * (1.0d0 - I_)      !Labor income
    elseif (gI_p_L==2) then  !leisure=1/2 part time
       rtemp = pw(git) * gH * (0.5d0 - I_)    !Labor income 
    endif
     
    call A_AIME_prime( c_, rtemp, rx2)
    if ((git < N_T) .and. (gcL_Iscapped)) then
       rx2(1) = min(sgridA%a(git_pos+1, VN_A(git+1)), rx2(1)) 
       rx2(2) = min(sgridAIME%a(git_pos+1, VN_AIME(git+1)), rx2(2))
    endif
    
    cmaxee_valfun_cI = - ( ufunc_c(c_) + cmaxee_ev((/I_,L_/), rx2) )
  end function cmaxee_valfun_cI

  !====================================================================
  ! Calculate value function given c, and I=gr_p_I, L=gr_p_L
  !====================================================================
  real(8) function cmaxee_valfun_c(rcx_)
    !calculate the value function given c for Period git
    ! to be minimized, so flip the sign
    REAL(8), INTENT(IN) :: rcx_(:)
    real(8) :: c_, I_, L_, rtemp, rx2(2)
    c_ = real2ab_exp(gr_p_cbound(1), gr_p_cbound(2), rcx_(1))
    I_ = gr_p_I
    L_ = gr_p_L    
    IF (c_ .NE. c_) THEN	! this will be true if I_ or L_ is NaN
       cmaxee_valfun_c = NR_BIG2
       return
    endif
    !The next period state variables
    if (gI_p_L==1) then !leisure = 1
       rtemp = 0.0d0
    elseif (gI_p_L==0) then  !leisure=0
       rtemp = pw(git) * gH * (1.0d0 - I_)      !Labor income
    elseif (gI_p_L==2) then  !leisure=1/2 part time
       rtemp = pw(git) * gH * (0.5d0 - I_)      !Labor income 
    endif
    
    call A_AIME_prime( c_, rtemp, rx2)
    if ((git < N_T) .and. (gcL_Iscapped)) then
        rx2(1) = min(sgridA%a(git_pos+1, VN_A(git+1)), rx2(1))
        rx2(2) = min(sgridAIME%a(git_pos+1, VN_AIME(git+1)), rx2(2))
    endif
    
    cmaxee_valfun_c = - ( ufunc_c(c_) + cmaxee_ev((/I_,L_/), rx2) )
  end function cmaxee_valfun_c
 
  !====================================================================
  ! Calculate value function given I, and c=grCF
  !====================================================================
  real(8) function cmaxee_valfun_I(rcx_)
    !calculate the value function given I for Period git
    ! to be minimized, so flip the sign
    REAL(8), INTENT(IN) :: rcx_(:)
    real(8) :: c_, I_, L_, rtemp, rx2(2)
    c_ = grCF
    I_ = real2ab_exp(gr_p_Ibound(gI_p_L,1), gr_p_Ibound(gI_p_L,2), rcx_(1))
    L_ = gr_p_L    
    IF (I_ .NE. I_) THEN	! this will be true if is NaN
       cmaxee_valfun_I = NR_BIG2
       return
    endif
    !The next period state variables
    if (gI_p_L==1) then !leisure = 1
       rtemp = 0.0d0
    elseif (gI_p_L==0) then  !leisure=0
       rtemp = pw(git) * gH * (1.0d0 - I_)      !Labor income
    elseif (gI_p_L==2) then  !leisure=1/2 part time
       rtemp = pw(git) * gH * (0.5d0 - I_)      !Labor income 
    endif
    
    call A_AIME_prime( 0.0d0, rtemp, rx2)
    rx2(1) = sgridAnrlb%a(git_pos+1) !A_prime in the case of c=grCF
    
    if ((git < N_T) .and. (gcL_Iscapped)) then
        rx2(1) = min(sgridA%a(git_pos+1, VN_A(git+1)), rx2(1))
        rx2(2) = min(sgridAIME%a(git_pos+1, VN_AIME(git+1)), rx2(2))
    endif
        
    cmaxee_valfun_I = - (ufunc_c(c_) + cmaxee_ev((/I_,L_/), rx2))
  end function cmaxee_valfun_I
  
  real(8) function cmaxee_ev(IL_, xp_)
    !calculate the continuation value, vnext (EMAX)
    REAL(8), INTENT(IN) :: IL_(:), xp_(:)  !xp_(1:2): A and AIME
    real(8) :: xnext(3), vnext(0:4), rtemp, rtemp2, rtemp3,  &
             & Aprime(0:N_FS_spinc_int), Hprime(0:N_XI_int)
    integer(4) :: ixnext(3), i, iFS_next, ihealth_next, jnextHS0, jnextHS1, j,  &
             & iAprime(0:N_FS_spinc_int), iHprime(0:N_XI_int)
        
    if (git >= N_T) then
       vnext = Vbequest_A(xp_(1))
    else
       !AIME prime 
       xnext(3) = xp_(2)
       ixnext(3) = locate(sgridAIME%a(git_pos+1,1:VN_AIME(git+1)), xnext(3))  !AIME prime
       if (ixnext(3) == 0) then
          ixnext(3) = 1
       elseif (ixnext(3) == VN_AIME(git+1) .AND. VN_AIME(git+1)>1) then
          ixnext(3) = VN_AIME(git+1) - 1
          if (gcL_Iscapped) xnext(3) = sgridAIME%a(git_pos+1,VN_AIME(git+1))
       endif

       !A prime
       do j = 0, N_FS_spinc_int
          if (j == 0) then
             Aprime(0) = xp_(1)  !A prime 
          else
             Aprime(j) = xp_(1) + gr_FS_spinc_xw(git+1,j,3) !A prime
          endif   
          iAprime(j) = locate(sgridA%a(git_pos+1,1:VN_A(git+1)), Aprime(j))  !A prime index
          if (iAprime(j) == 0) then
             iAprime(j) = 1
          elseif (iAprime(j) == VN_A(git+1) .AND. VN_A(git+1)>1) then
             iAprime(j) = VN_A(git+1) - 1
             if (gcL_Iscapped) Aprime(j) = sgridA%a(git_pos+1,VN_A(git+1))
          endif
       enddo !do j

       !H prime
       if (IL_(2)>=0.95d0 .AND. (gcI_robust .NE. 8 .AND. gcI_robust .NE. 10 .AND. gcI_robust .NE. 12)) then  !if not working then not hit by the shock, except exo case
          j = 0 
       else
          j = N_XI_int
       endif
       do i = 0, j
          if (i == 0) then  !not working
             Hprime(i) = (1.0d0 - psigma_offjob)*gH   !H prime 
          else  !if working or exo
             if (IL_(2)>=0.95d0) then 
                 rtemp = Hprime(0)  !not working (exo cases only)
             else
                 rtemp = (1.0d0 - psigma_onjob)*gH   !H prime, on job
             endif
             
             if (gcI_robust == 8 .OR. gcI_robust == 9) then
                Hprime(i) = rtemp + max(0.0d0,  &
                    &  gr_XI_xw(i,3)*ppi*(1.0d0+palpha_I*dble(git-N_0)/gr_LBD_t+palpha_H*dble((git-N_0)**2)/gr_LBD_tsq))
             elseif (gcI_robust == 10 .OR. gcI_robust == 11) then  !exo vH (regardless of working) lbd vH (working)
                Hprime(i) = rtemp + max(0.0d0,  &
                    &  gr_XI_xw(i,3)*ppi*(1.0d0 + palpha_I*gH + palpha_H*gH**2 ))
             elseif (gcI_robust == 12 .OR. gcI_robust == 13) then  !exo Hone (regardless of working) lbd Hone (working)
                Hprime(i) = rtemp + max(0.0d0, gr_XI_xw(i,3)*ppi*(gH**palpha_H))    
             else
                Hprime(i) = rtemp + gr_XI_xw(i,3)*ppi*(IL_(1)**palpha_I)*(gH**palpha_H)
             endif
          endif    !if (i == 0)
          iHprime(i) = locate(sgridH%a(git_pos+1,:), Hprime(i))  !H prime index
          if (iHprime(i) == 0) then
             iHprime(i) = 1
          elseif (iHprime(i) == N_H) then
             iHprime(i) = N_H - 1
             if (gcL_Iscapped) Hprime(i) = sgridH%a(git_pos+1,N_H)
          endif
       enddo !do i
       
       vnext = 0.0d0
       !health status next period, starting from IHasHealthAge-1.
       !Starting from IHasHealthAge-1, inclusive, for the NEXT period (git+1) there will be multiple health status
       !in the baseline model without health, or models where health starts from age IHasHealthAge
       if (IHasHealthAge == 0) then   !no health
          jnextHS0 = 1
          jnextHS1 = 1
       elseif (git<(IHasHealthAge-1)) then
          jnextHS0 = 2
          jnextHS1 = 2
       elseif (N_HEALTH>1 .AND. giHEALTH==N_HEALTH) then !DI is absorbing
          jnextHS0 = N_HEALTH
          jnextHS1 = N_HEALTH
       else
          jnextHS0 = 1
          jnextHS1 = N_HEALTH
       endif       
       do ihealth_next = jnextHS0, jnextHS1  !health status next period
          if (IL_(2)>=0.95d0 .AND. (gcI_robust .NE. 8 .AND. gcI_robust .NE. 10 .AND. gcI_robust .NE. 12)) then  !if not working then not hit by the shock
             !H prime
             xnext(2) = Hprime(0) 
             ixnext(2) = iHprime(0)
             !Family status next period, single (1) or married spouse not working (2)
             xnext(1) = Aprime(0)  !A prime
             ixnext(1) = iAprime(0)
             do iFS_next = 1, N_FS-1 !single (1) or married spouse not working (2)
                if (gcL_triangle) then 
                   rtemp = linear_inte_tetrahedra3(sgridA%a(git_pos+1,ixnext(1):(ixnext(1)+1)), &
                      &  sgridAIME%a(git_pos+1,ixnext(3):(ixnext(3)+1)),sgridH%a(git_pos+1,ixnext(2):(ixnext(2)+1)), &
                      &  dvV_EV%a(ixnext(1):(ixnext(1)+1),ixnext(3):(ixnext(3)+1),ihealth_next, &
                      &  giRECSS_next,iFS_next,ixnext(2):(ixnext(2)+1),git_pos+1),  &
                      &  (/xnext(1),xnext(3),xnext(2)/))  !EV 
                else
                   rtemp = linear_inte3(gcI_inte_v_log, (/1,1,1/), sgridA%a(git_pos+1,ixnext(1):(ixnext(1)+1)), &
                      &  sgridAIME%a(git_pos+1,ixnext(3):(ixnext(3)+1)),sgridH%a(git_pos+1,ixnext(2):(ixnext(2)+1)), &
                      &  dvV_EV%a(ixnext(1):(ixnext(1)+1),ixnext(3):(ixnext(3)+1),ihealth_next, &
                      &  giRECSS_next,iFS_next,ixnext(2):(ixnext(2)+1),git_pos+1),  &
                      &  (/xnext(1),xnext(3),xnext(2)/))  !EV
                endif   
                vnext(ihealth_next) = vnext(ihealth_next) + gr_FS_trans(git,giFS,iFS_next)*rtemp
             enddo !do iFS_next
             !Family status next period, married and spouse working (3)
             !integrate over the spouse income next period
             iFS_next = 3
             rtemp2 = 0.0d0
             do j = 1, N_FS_spinc_int
                xnext(1) = Aprime(j) !A prime
                ixnext(1) = iAprime(j)  !A prime index
                if (gcL_triangle) then 
                   rtemp = linear_inte_tetrahedra3(sgridA%a(git_pos+1,ixnext(1):(ixnext(1)+1)), &
                      &  sgridAIME%a(git_pos+1,ixnext(3):(ixnext(3)+1)),sgridH%a(git_pos+1,ixnext(2):(ixnext(2)+1)), &
                      &  dvV_EV%a(ixnext(1):(ixnext(1)+1),ixnext(3):(ixnext(3)+1),ihealth_next,giRECSS_next, &
                      &  iFS_next,ixnext(2):(ixnext(2)+1),git_pos+1),  &
                      &  (/xnext(1),xnext(3),xnext(2)/))  !EV 
                else
                   rtemp = linear_inte3(gcI_inte_v_log, (/1,1,1/), sgridA%a(git_pos+1,ixnext(1):(ixnext(1)+1)), &
                      &  sgridAIME%a(git_pos+1,ixnext(3):(ixnext(3)+1)),sgridH%a(git_pos+1,ixnext(2):(ixnext(2)+1)), &
                      &  dvV_EV%a(ixnext(1):(ixnext(1)+1),ixnext(3):(ixnext(3)+1),ihealth_next,giRECSS_next, &
                      &  iFS_next,ixnext(2):(ixnext(2)+1),git_pos+1),  &
                      &  (/xnext(1),xnext(3),xnext(2)/))  !EV
                endif   
                rtemp2 = rtemp2 + gr_FS_spinc_xw(git+1,j,2)*rtemp
             enddo !do j
             vnext(ihealth_next) = vnext(ihealth_next) + gr_FS_trans(git,giFS,iFS_next)*rtemp2 
          else  !if working, then hit by the shock; or gcI_robust==8
             !! Gauss-Hermite integration for H
             do i=1, N_XI_int
                xnext(2) = Hprime(i)
                ixnext(2) = iHprime(i)
                !Family status next period, single (1) or married spouse not working (2)
                rtemp3 = 0.0d0
                xnext(1) = Aprime(0)  !A prime
                ixnext(1) = iAprime(0)
                do iFS_next = 1, N_FS-1
                   if (gcL_triangle) then 
                      rtemp = linear_inte_tetrahedra3(sgridA%a(git_pos+1,ixnext(1):(ixnext(1)+1)), &
                         &  sgridAIME%a(git_pos+1,ixnext(3):(ixnext(3)+1)),sgridH%a(git_pos+1,ixnext(2):(ixnext(2)+1)), &
                         &  dvV_EV%a(ixnext(1):(ixnext(1)+1),ixnext(3):(ixnext(3)+1),ihealth_next,giRECSS_next, &
                         &  iFS_next,ixnext(2):(ixnext(2)+1),git_pos+1),  &
                         &  (/xnext(1),xnext(3),xnext(2)/))  !EV 
                   else
                      rtemp = linear_inte3(gcI_inte_v_log, (/1,1,1/), sgridA%a(git_pos+1,ixnext(1):(ixnext(1)+1)), &
                         &  sgridAIME%a(git_pos+1,ixnext(3):(ixnext(3)+1)),sgridH%a(git_pos+1,ixnext(2):(ixnext(2)+1)), &
                         &  dvV_EV%a(ixnext(1):(ixnext(1)+1),ixnext(3):(ixnext(3)+1),ihealth_next,giRECSS_next, &
                         &  iFS_next,ixnext(2):(ixnext(2)+1),git_pos+1),  &
                         &  (/xnext(1),xnext(3),xnext(2)/))  !EV
                   endif   
                   rtemp3 = rtemp3 + gr_FS_trans(git,giFS,iFS_next)*rtemp
                enddo !do iFS_next
                !Family status next period, married and spouse working (3)
                !integrate over the spouse income next period
                iFS_next = 3
                rtemp2 = 0.0d0             
                do j = 1, N_FS_spinc_int
                   xnext(1) = Aprime(j)  !A prime
                   ixnext(1) = iAprime(j)  !A prime index
                   if (gcL_triangle) then 
                      rtemp = linear_inte_tetrahedra3(sgridA%a(git_pos+1,ixnext(1):(ixnext(1)+1)), &
                         &  sgridAIME%a(git_pos+1,ixnext(3):(ixnext(3)+1)),sgridH%a(git_pos+1,ixnext(2):(ixnext(2)+1)), &
                         &  dvV_EV%a(ixnext(1):(ixnext(1)+1),ixnext(3):(ixnext(3)+1),ihealth_next,giRECSS_next, &
                         &  iFS_next,ixnext(2):(ixnext(2)+1),git_pos+1),  &
                         &  (/xnext(1),xnext(3),xnext(2)/))  !EV 
                   else
                      rtemp = linear_inte3(gcI_inte_v_log, (/1,1,1/), sgridA%a(git_pos+1,ixnext(1):(ixnext(1)+1)), &
                         &  sgridAIME%a(git_pos+1,ixnext(3):(ixnext(3)+1)),sgridH%a(git_pos+1,ixnext(2):(ixnext(2)+1)), &
                         &  dvV_EV%a(ixnext(1):(ixnext(1)+1),ixnext(3):(ixnext(3)+1),ihealth_next,giRECSS_next, &
                         &  iFS_next,ixnext(2):(ixnext(2)+1),git_pos+1),  &
                         &  (/xnext(1),xnext(3),xnext(2)/))  !EV
                   endif   
                   rtemp2 = rtemp2 + gr_FS_spinc_xw(git+1,j,2)*rtemp
                enddo !do j
                rtemp3 = rtemp3 + gr_FS_trans(git,giFS,iFS_next)*rtemp2
                vnext(ihealth_next) = vnext(ihealth_next) + gr_XI_xw(i,2) * rtemp3
             enddo !do i=1, N_XI_int
          endif !if (IL_(2)>=0.95d0 .AND. (gcI_robust .NE. 8))
       enddo !do ihealth_next
    endif  !if (git >= N_T)
    
    !health transition probability gr_health_matrix(N_0:N_T,1:4,0:4)
    !the loop control makes sure that for IHasHealthAge>1, it is always git>=IHasHealthAge for giHEALTH==2
    if (git >= N_T) then
       cmaxee_ev = pbeta * vnext(1)
    elseif (IHasHealthAge == 0) then   !no health
       cmaxee_ev = pbeta * vnext(1)
    elseif (git<(IHasHealthAge-1)) then !has health but before "IHasHealthAge-1"
       cmaxee_ev = pbeta * vnext(2)
    elseif (N_HEALTH>1 .AND. giHEALTH==N_HEALTH) then !DI is absorbing
       cmaxee_ev = pbeta * vnext(giHEALTH) 
    else
       if ( (gcI_CF==41 .OR. gcI_CF==46) .AND. git>=50 ) then 
          !41.	Before 50 same as the baseline model of health, then the health status at age 50 is fixed for the rest of the life-cycle
          !46.	Before 50 same as the baseline model of health, then the health status at age 50 is fixed for the rest of the life-cycle, 
          !     and the taste for leisure stays same as age 50 and does not increase with the age.
          cmaxee_ev = pbeta * vnext(giHEALTH) 
       elseif ( (gcI_CF==42 .OR. gcI_CF==47) .AND. git>=49) then 
          !42.	Before 50 same as the baseline model of health, then the health status at age 50 becomes excellent for everyone, 
          !     and then stays excellent for the rest of the life-cycle. 
          !47.	Before 50 same as the baseline model of health, then the health status at age 50 becomes excellent for everyone, 
          !     and then stays excellent for the rest of the life-cycle,
          !     and the taste for leisure stays same as age 50 and does not increase with the age.  
          cmaxee_ev = pbeta * vnext(1)
       elseif ( (gcI_CF==43 .OR. gcI_CF==48) .AND. git>=49) then 
          !43.	Before 50 same as the baseline model of health, then the health status at age 50 becomes good for everyone, 
          !     and then stays good for the rest of the life-cycle. 
          !48.	Before 50 same as the baseline model of health, then the health status at age 50 becomes good for everyone, 
          !     and then stays good for the rest of the life-cycle,
          !     and the taste for leisure stays same as age 50 and does not increase with the age.  
          cmaxee_ev = pbeta * vnext(2)   
       elseif ((gcI_CF==44 .OR. gcI_CF==49) .AND. git>=49) then 
          !44.	Before 50 same as the baseline model of health, then the health status at age 50 becomes bad for everyone, 
          !     and then stays bad for the rest of the life-cycle. 
          !49.	Before 50 same as the baseline model of health, then the health status at age 50 becomes bad for everyone, 
          !     and then stays bad for the rest of the life-cycle, and the taste for leisure stays same as age 50 and does not increase with the age. 
          cmaxee_ev = pbeta * vnext(3)
       elseif ((gcI_CF==45 .OR. gcI_CF==50) .AND. git>=49) then 
          !45.	Before 50 same as the baseline model of health, then the health status at age 50 becomes disabled for everyone, 
          !     and then stays disabled for the rest of the life-cycle. 
          !50.	Before 50 same as the baseline model of health, then the health status at age 50 becomes disabled for everyone, 
          !     and then stays disabled for the rest of the life-cycle, and the taste for leisure stays same as age 50 and does not increase with the age. 
          cmaxee_ev = pbeta * vnext(4)   
       else  !the model with health status
          vnext(0) = 0.0d0
          do j = 1, 4
             if (IHasHealthAge>0 .AND. git==(IHasHealthAge-1)) then
                rtemp = gr_health_matrix(git+1,j,0)  !prob of having j health status next period
             else
                rtemp = gr_health_matrix(git,giHEALTH,j)
             endif  !if (IHasHealthAge>0)
             vnext(0) = vnext(0) + rtemp*vnext(j)
          enddo !do j = 1, 4
          cmaxee_ev = pbeta * vnext(0)
       endif      !if (IHasHealthAge == 0) 
    endif   !if (git >= N_T)
    
    return
  END function cmaxee_ev

  !====================================================================
  ! Calculate the value given I and L, for the last period (N_T)
  ! I_ = 0
  ! directly using EE
  ! consumption is stored in gr_cc
  !====================================================================
  real(8) function cmaxee_valfun_pT(iL_)
    integer, INTENT(IN) :: iL_
    real(8) :: I_, Vrcx2(2), rtemp, rL
        
    IF (iL_ .NE. iL_) THEN	! this will be true if I_ or L_ is NaN
       cmaxee_valfun_pT = - NR_BIG2
       return
    endif
    
    I_ = 0.0d0
      
    
    !The next period state variables. No need for H since this is the last period.
    if (iL_ == 1) then !leisure==1
       rtemp = 0.0d0
       rL = 1.0d0
    elseif (iL_ == 0) then
       rtemp = pw(git) * gH      !Labor income
       rL = 0.0d0
    else  !part time, iL_=2
       rtemp = pw(git) * gH * 0.5d0
       rL = 0.5d0
    endif
    
    call A_AIME_prime( 0.0d0, rtemp, Vrcx2) !(c,wageinc,/Aprime,AIMEprime/)
    
    
    ! solve optimal consumption, for the last period
    !Vrcx2(1) is Aprime for I, L and c=0, or Aprime_c0
    gr_cc = gr_bequest_AP(giFS,iL_) * (pb2 + Vrcx2(1))    
    
    !the asset next period
    rtemp = Vrcx2(1) - gr_cc  !A_prime
    if (rtemp < sgridAnrlb%a(N_T_grid+1)) then
       gr_cc = Vrcx2(1) - sgridAnrlb%a(N_T_grid+1)
       Vrcx2(1) = sgridAnrlb%a(N_T_grid+1) !A_prime
    else
       Vrcx2(1) = rtemp  !A_prime
    endif
    
    if (gr_cc > 0.0d0) then
       cmaxee_valfun_pT = ufunc_c(gr_cc) + cmaxee_ev((/I_,rL/), Vrcx2)
    else
       cmaxee_valfun_pT = -NR_BIG2 
    endif
    
    !check c=grCF case
    if ( min(gA,gr_cc) < 20.0d0 ) then
       !Vrcx2(1) = sgridA%a(git,1)*(1.0d0+pr)  !here git==N_T
       Vrcx2(1) = sgridAnrlb%a(N_T_grid+1)
       rtemp = ufunc_c(grCF) + cmaxee_ev((/I_,rL/), Vrcx2)
       if (rtemp > cmaxee_valfun_pT) then
          gr_cc = grCF
          cmaxee_valfun_pT = rtemp
       endif
    endif !if ( min(gA,gr_cc) < 20.0d0 ) then 
       
    return
  END function cmaxee_valfun_pT

  
  
  !====================================================================
  !===== A prime and AIME prime
  ! It is better not to be separated into two function because in calculating Aprime
  ! we have to take into account the earnings test
  ! given git, gA, gAIME, giRECSS, giRECSS_next (current policy) 
  subroutine A_AIME_prime( c_, winc_, aaprime_ )
     IMPLICIT NONE
     real(8), intent(in) :: c_, winc_  !consumption, wage income
     real(8), intent(out) :: aaprime_(:)  !A prime, AIME prime
     integer :: diffage 
     real(8) :: ssinc, ssdiinc, pia, piaprime, rage, aimeadj_, rothinc, rtemp, rtemp2, rA
     
     aimeadj_ = 1.0d0     
     ssinc = 0.0d0
     ssdiinc = 0.0d0
     pia = AIME2PIA(gAIME)  !Effective PIA and  effective AIME     

     !===========================================================================================================
     ! Other income (on age, agesq, dage62, dage65, col, constant)
     rothinc = 0.0d0
          
     !====================================================================
     ! Social security rules:
     ! 1. early retirement
     ! 2. delayed retirement credit 
     ! output: ssinc, 
     !         current_AIMEprime_adj (use this to adjust AIME next period)
     !=====================================================================
     
     !===========================================================================
     ! In most cases, that��s going to be disability. Regardless of your age 
     ! when you start getting Social Security Disability Insurance (SSDI), 
     ! you receive what you would get if you claimed benefits at full retirement age (FRA) 
     ! �� the age at which you are entitled to 100 percent of the benefit calculated 
     ! from your average monthly earnings. (FRA Is currently 66 and rising gradually 
     ! to 67 for people born in 1960 or later.) 
     ! If you turn 62 in 2019, you��re eligible for only 72.5 percent of that retirement 
     ! benefit, so you��re probably better off sticking with SSDI. When you reach FRA, 
     ! the disability benefit automatically converts to a retirement benefit, and you��ll 
     ! get the same monthly amount you��ve been getting.
     !===========================================================================
     if (giHEALTH==4 .AND. git<pNRA .AND. winc_<gr_SSDI_SGA) then  !SSDI and SSI
        ssdiinc = pia * 12.0d0   !translate monthly PIA to yearly SSDI income
        ssdiinc = max(ssdiinc, gr_SSI)
        if ((git - N_0) >= 35) then
            rtemp = 1.0d0
        else
            rtemp = 0.0d0
        endif
        
        !update AIME prime
        if (gcI_CF == 26) then      !26---Modifying C5: AIME_t+1 = w_t (last year only) if L_t<1
            if (winc_ <= 1.0d-8) then
                aaprime_(2) = gAIME
            else
                aaprime_(2) = min(winc_,gr_sse_cap)/12.0d0
            endif
        elseif (gcI_CF == 27) then  !27---Modifying C5: AIME_t+1 = max(AIME_t, w_t) (highest year only)
            aaprime_(2) = max(gAIME, min(winc_,gr_sse_cap)/12.0d0)
        else
            aaprime_(2) = gAIME + max(0.0d0, min(winc_,gr_sse_cap)/420.0d0 - rtemp*cps_aimeshare%a(git)*gAIME)
        endif    
     elseif (giRECSS_next == 1) then  ! do not receive ss yet in current (and next) period (giRECSS_next == 1)
        ssinc = 0.0d0
        !AIME prime 35*12=420. winc_ is the yearly wage income here.
        !update AIME
        if ((git - N_0) >= 35) then
            rtemp = 1.0d0
        else
            rtemp = 0.0d0
        endif
        !update AIME prime
        if (gcI_CF == 26) then      !26---Modifying C5: AIME_t+1 = w_t (last year only) if L_t<1
            if (winc_ <= 1.0d-8) then
                aaprime_(2) = gAIME
            else
                aaprime_(2) = min(winc_,gr_sse_cap)/12.0d0
            endif
        elseif (gcI_CF == 27) then  !27---Modifying C5: AIME_t+1 = max(AIME_t, w_t) (highest year only)
            aaprime_(2) = max(gAIME, min(winc_,gr_sse_cap)/12.0d0)
        else
            aaprime_(2) = gAIME + max(0.0d0, min(winc_,gr_sse_cap)/420.0d0 - rtemp*cps_aimeshare%a(git)*gAIME)
        endif    
     elseif (giRECSS_next == 2) then  !receiving ss
        if (giRECSS == 2) then  !(giRECSS == 2) & (giRECSS_next == 2)
           ! has already taken ss in previous period
           ssinc = pia * 12.0d0 !translate monthly PIA to yearly ss income 
        elseif (giRECSS == 1) then
           ! didn't take ss in previous period
           ! current period is the first time to take ss
           if (git < pERA) then
              ssinc = 0.0d0
           elseif (git < pNRA) then  !early retirement
              diffage = pNRA - git
              if (diffage <= 3) then
                 piaprime = pia * (1.0d0 - 0.0667d0 * dble(diffage))
              else
                 piaprime = pia * (0.8d0 - 0.05d0 * dble(diffage - 3))
              endif
              ssinc = piaprime * 12.0d0 !translate monthly PIA to yearly ss income
              rtemp = PIA2AIME( piaprime )  !new AIME
              if (gAIME > 0.0d0) aimeadj_ = rtemp / gAIME
           elseif ( gcL_IsDRC .AND. (git > pNRA) ) then  !DRC: delayed retirement credit
              rtemp = 1.0d0+dble(max(0, min(git,69)-pNRA))*0.06d0
              piaprime = pia * rtemp
              ssinc = piaprime * 12.0d0 !translate monthly PIA to yearly ss income
              rtemp = PIA2AIME( piaprime )
              if (gAIME > 0.0d0) aimeadj_ = rtemp / gAIME
           else  !retire at NRA or (after NRA and no DRC)
              ssinc = pia * 12.0d0 !translate monthly PIA to yearly ss income 
           endif  !if (git < pERA)
        endif    !(giRECSS == 2)
        
        !====================================================================
        ! still (giRECSS_next == 2) receiving SS
        ! earnings test (including Delayed retirement credit from withdrew ET)
        !!1 ---Removing Social Security earnings test.
        !=====================================================================
        if ((gcI_CF .NE. 1) .AND. (git>=pERA) .AND. (git<70)) then !there is earnings test
           if ( ((git<pNRA) .and. (winc_>=pETESTbp(1))) .or. ((git>=pNRA) .and. (winc_>=pETESTbp(2))) ) then  !before 2000
              if (git<pNRA) then
                 rtemp = winc_ - pETESTbp(1) 
                 rtemp = min(ssinc, rtemp / 2.0d0)  !withheld benefit
              else                 
                 rtemp = winc_ - pETESTbp(2) 
                 rtemp = min(ssinc, rtemp / 3.0d0)  !withheld benefit
              endif
              if (gcL_IsDRC .and. (rtemp > 0.0d0) .and. (ssinc > 0.0d0)) then  !!Delayed retirement credit
                 piaprime = gAIME * aimeadj_  !AIME adjusted
                 piaprime = AIME2PIA(piaprime) !PIA prime
                 if (git<pNRA) then
                    piaprime = piaprime * (1.0d0 + 0.067d0 * (rtemp / ssinc) )
                 else
                    piaprime = piaprime * (1.0d0 + 0.06d0 * (rtemp / ssinc) ) 
                 endif
                 piaprime = PIA2AIME(piaprime) !AIME prime
                 if (gAIME > 0.0d0) aimeadj_ = piaprime / gAIME
              endif !if (gcL_IsDRC .and. (rtemp > 0.0d0))
              ssinc = max(0.0d0, ssinc - rtemp)
           endif !if ( ((git<pNRA)
        endif !if (gcI_CF .NE. 1
        
        aaprime_(2) = gAIME * aimeadj_  !no update from wage income after starting claiming SSB
     endif ! if (giHEALTH==4 .AND. git<pNRA .AND. winc_<gr_SSDI_SGA)

     !3.No Social Security system (no SS tax, no SS benefit).
     !6.Remove SS benefit, but keep SS tax.
     !5.Reduce SS benefit by 20%. SS tax unchanged.
     !10.No Social Security system (no SS tax, no SS benefit)---keep SSDI unchanged
     !11.Reduce SS benefit by 20%. SS tax unchanged.---keep SSDI unchanged
     !12.Remove SS benefit, but keep SS tax;---keep SSDI unchanged
     if ( (gcI_CF==3) .OR. (gcI_CF==6) ) then
        ssinc = 0.0d0
        ssdiinc = 0.0d0
     elseif (gcI_CF == 5) then
        ssinc = 0.8d0 * ssinc
        ssdiinc = 0.8d0 * ssdiinc
     elseif ( (gcI_CF==10) .OR. (gcI_CF==12) ) then
        ssinc = 0.0d0
     elseif (gcI_CF == 11) then
        ssinc = 0.8d0 * ssinc   
     endif     
     
     !=================================
     ! A prime
     !=================================
     !aaprime_(1) = (1.0d0+pr)*gA + ssinc + winc_ + rothinc - c_  ! A prime
     !aaprime_(1) = winc_ + rothinc  !capital interest income is NOT taxable
     !capital interest income is taxable
     rA = pr*gA
     aaprime_(1) = winc_ + rothinc + max(0.0d0, rA)  
     !taxable social security benefit
     if (aaprime_(1) <= 0.0d0) then
        aaprime_(1) = ssinc
     else
        rtemp = aaprime_(1) + (ssinc / 2.0d0)
        if (rtemp <= gr_sstax_1) then
           rtemp2 = 0.0d0 !taxable ss benefit 
        else    
           rtemp2 = min(0.85d0*ssinc, 0.5d0*min(ssinc, rtemp-gr_sstax_1, gr_sstax_2)+0.85d0*max(0.0d0, rtemp-gr_sstax_3))  !taxable ss benefit
        endif   
        rtemp = ssinc - rtemp2  !nontaxable ss part
        aaprime_(1) = ftaxcode(aaprime_(1)+rtemp2) + rtemp  !after tax income
     endif   
     
     !A prime, positive capital interest income is taxable
     aaprime_(1) = gA + min(0.0d0, rA) + aaprime_(1) - c_
     aaprime_(1) = aaprime_(1) + ssdiinc  !plus SSDI & SSI
     
  end subroutine A_AIME_prime


  !========================================================================================================
  ! Check consistency for V and C, for given  (:,giAIME,giHEALTH,giRECSS,giFS,giH,git), Along A
  ! We are not using Euler Equation to solve optimal C, but we will interpolate C in simulation
  !========================================================================================================
  SUBROUTINE cmaxee_consistency_VCinA()
     IMPLICIT NONE
     integer(4), parameter :: j = 3  !only the 3 grids in the lower end.
     integer(4) :: i
     real(8) :: Vrtemp(j), rtemp
     
     !It is important that at the lower end, V is strictly increasing with A, to avoid the extrapolation problem.     
     do i = 0, Ipt+1
        !V 
        Vrtemp(1:j) = dvV(i)%a(1:j,giAIME,giHEALTH,giRECSS,giFS,giH,git_pos)
        if (Vrtemp(1) >= Vrtemp(2)) then !need to fix it
           call sort(Vrtemp(1:j)) 
           dvV(i)%a(1:j,giAIME,giHEALTH,giRECSS,giFS,giH,git_pos) =  &
              &  dvV(i)%a(j,giAIME,giHEALTH,giRECSS,giFS,giH,git_pos)-Vrtemp(j)+Vrtemp(1:j) 
        endif !if (Vrtemp(1) >= Vrtemp(2))
        
     enddo ! do i     
  end subroutine cmaxee_consistency_VCinA
   
  !========================================================================================================
  ! Check consistency for EV, for given  (:,giAIME,giHEALTH,giRECSS,giFS,giH,git), Along A
  !========================================================================================================
  SUBROUTINE cmaxee_consistency_EVinA()
     IMPLICIT NONE
     integer(4), parameter :: j = 3  !only the 3 grids in the lower end.
     real(8) :: Vrtemp(5), rtemp
     
     !It is important that at the lower end, V is strictly increasing with A, to avoid the extrapolation problem.     
     Vrtemp(1:j) = dvV_EV%a(1:j,giAIME,giHEALTH,giRECSS,giFS,giH,git_pos)
     if (Vrtemp(1) >= Vrtemp(2)) then !need to fix it
        call sort(Vrtemp(1:j)) 
        dvV_EV%a(1:j,giAIME,giHEALTH,giRECSS,giFS,giH,git_pos) =  &
            &   dvV_EV%a(j,giAIME,giHEALTH,giRECSS,giFS,giH,git_pos)-Vrtemp(j)+Vrtemp(1:j)
     endif !if (Vrtemp(1) >= Vrtemp(2))
  end subroutine cmaxee_consistency_EVinA
  
  
END MODULE CMAXEE
